Critical Temperature and Thermodynamics of Attractive Fermions at Unitarity 
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The unitarity regime of the BCS-BEC crossover can be realized by diluting a system of two- 
component lattice fermions with an on-site attractive interaction. We perform a systematic-error- 
free finite-temperature simulations of this system by diagrammatic determinant Monte Carlo. The 
critical temperature in units of Fermi energy is found to be T c /sf = 0.152(7). We also report the 
behavior of the thermodynamic functions, and discuss the issues of thermometry of ultracold Fermi 
gases. 
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PACS numbers: 03.75.Ss,05.10.Ln,71.10.Fd 

The unitarity limit is commonly referred to as the limit 
of a diverging scattering length a — > oo, and an effective 
range of the interaction r e — > 0. A Fermi gas in this 
limit attains universality: at low enough temperature the 
only relevant length scale is given by the density, n, since 
the divergent scattering length drops out completely and 
the system's properties are independent of the interac- 
tion details. The unitarity limit is approximately real- 
ized in the inner crust of the neutron stars, where the 
neutron-neutron scattering length is nearly an order of 
magnitude larger than the mean interparticle separation 
Unitarity conditions can also be achieved with cold 
trapped atom gases using the Feshbach resonance tech- 
nique, i.e. tuning the scattering length to infinity using 
the magnetic field. In recent years these systems have 
been extensively studied experimentally 0, 0] • 

In the limit of £ = 1/ na — > +oo the fermions pair into 
bosonic molecules and form a Bose-Einstein condensate 
(BEC). In the opposite limit £ — > — oo one recovers the 
Bardeen-Cooper-Schrieffer (BCS) limit. The unitarity 
limit £ — > separates these two extremes. In all these 
cases, a gas undergoes a superfluid (SF) phase transition 
at some temperature, which depends on £. 

The early analytical treatments of the unitary Fermi 
gas have been based on the extension of the BCS-type 
many-body wave function Most of the subsequent 
elaborations are also of mean-field type (with or with- 
out fluctuations) @, @, 0, @, S E3] • The accuracy and 
reliability of such approximations is nevertheless ques- 
tionable given the strongly interacting nature of the uni- 
tarity regime, and the results differ by nearly an order of 
magnitude. 

Monte Carlo (MC) simulations of Fermi systems are, 
in general, severely hindered by a sign problem ll| . For- 
tunately for fermions with attractive contact interaction 
the sign problem can be avoided 12|, [l3j, [l4[ . The ground 
state of a unitary Fermi gas has been studied within a 
fixed- node MC framework [HI, the systematic errors of 
which depend on the quality of a guess of the nodal struc- 
ture of a many-body fermion configuration. Despite a 



number of calculations at finite temperatures 1(| 17, 18| . 
a reliable estimate of the critical temperature is lack- 
ing. The purpose of this Letter is to provide accurate 
results for the critical temperature and thermodynamic 
functions of a three-dimensional (3D) unitary Fermi gas 
using a novel determinant diagrammatic MC method free 
of systematic errors. 

Consider an attractive Hubbard model (AHM) defined 
by the Hamiltonian H = Hq + H\ , with 



Hn = 



ICCT 



Cko 



Hi = -U 2^n xT n xi , (1) 



= A 



where c kl7 is a fermion creation operator, 
a =t, I is the spin index, x enumerates sites of a 3D 
simple cubic lattice with periodic boundary conditions, 
the quasimomentum k spans the corresponding Brillouin 
zone, ek = — 2iJ^_ 1 cosfc a is the tight-binding spec- 
trum, t = 1 is the the hopping amplitude, (i stands for 
the chemical potential, U > is the on-site attraction, 
and we have set the lattice spacing to unity. 

By solving the two-body problem of the model (fT]) 
one finds that the scattering length diverges at U c — 
(L- 3 EkeBZ V^feT 1 ~ 7.915*. We use this value of U 
throughout. Since we are ultimately interested in the 
continuum rather than lattice results, we study the low 
density limit v — > 0, where ^ v ^ 2 is the filling frac- 
tion. We define Fermi momentum as kp = (St: 2 ^) 1 ^ 3 and 
Fermi energy ep = k F , as those of a continuum gas with 
the same effective mass and number density n = v. 

We simulate the model (UJ by diagrammatic deter- 
minant MC, discussed in detail in Refs. [l3|, One 
starts by expanding cxp(— (3H) in the interaction repre- 
sentation in powers of Hi . The resulting Feynmann dia- 
grams consist of four-point vertices representing the Hub- 
bard interaction, connected by free single-particle prop- 
agators. The sum over all possible ways of connecting 
vertices with propagators, in the n-th order diagram is 
represented by a vertex configuration S n — {(x^, Tj),j = 
1, ...,n)}, where r is the imaginary time, see Fig. [TJ . 
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FIG. 1: A sketch of a vertex configuration for the correla- 
tion function. Brown dots are the four-point vertices, with 
the incoming and outgoing lines shown. Red diamonds rep- 
resent the two-point vertices corresponding to P(x, r) and 
P'(x',t'). See the text for discussion. 
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In case of equal number of spin-up and spin-down parti- 
cles, the differential weight of a configuration is positive 
definite: 



dV(S n ) = U n \ det A(5„)| 2 Yl d Tj , 



(2) 



where A(iS„) is an n x n matrix built on single-particle 
propagators: Aij — G^(^i — Xj,Ti —Tj). 

The configuration space is sampled with worm-type 
pol ] updating scheme [2l[ , based on the two-particle cor- 
relation function 



G 2 (x,r;x',r') = (T r F(x,r)pt( x ', T '))^ 



(3) 



where T T is the r— ordering, P(x, r) = c x ^(t)c x i(t) is the 
pair annihilation operator, a normalization factor J\f — 
f3L 3 is introduced for future convenience (/3 is an inverse 
temperature), and (■ ■ • } is the thermal average. The non- 
zero asymptotic value of// drdr / G , 2(x, r; x', r') as |x — 
x' | — > oo is proportional to the condensate density. 

Fig. [1] shows a sketch of a vertex configuration, which 
features a pair of two-point vertices associated with 
P(x,t) and P(x',r'). 

The typical number of vertices in a configuration 
scales with the system volume as M oc (3UL 3 . Thus 
the Metropolis acceptance ratios for the updates in- 
volve the ratio of macroscopically large determinants 
det A(S' n , )/det A(S n ) with n' = n or n ± 1. Since we 
only need ratios of determinants, fast-update formulas 
fl3f can be used to reduce the computational complexity 
of an update from M 3 down to M 2 . 

We validate our method by comparing results against 
the exact diagonalization data for a 4 x 4 cluster [22j ]. 
and simulations of a critical temperature at quarter filling 
H [13. In both cases we find agreement within a few 
percent accuracy. 

We work in the grand canonical ensemble at fixed 
(L, T. fi). Extracting the unitarity limit critical tem- 
perature, T c , for the continuum gas from the lattice sim- 



FIG. 2: A typical crossing of the R(L, T) curves. The er- 
rorbars are 2a, and solid lines are the linear fits to the MC 
points. The data are for the \i = —5.2t, thus v(T = T C ,L — > 
oo) = 0.148(1). 



ulation is a two-stage process: first, we study the ther- 
modynamic limit L — y oo to obtain T c (y) at a given v 
,and then extrapolate to the continuum limit v — y 0. 

The first task is performed as follows: we simu- 
late a series of system sizes L\ > L% > . . . at vari- 
ous temperatures. At the critical point the correlation 
function ([3]) decays as a power-law at large distances: 
G2(x, r; x', r') cx l/|x — x'| 1+? ', where 77 is an anomalous 
dimension. Since the transition is expected to belong to 
the U(l) universality class, we use 77 = 0.038 |25| . Hence, 
if one sums and rescales the correlation function ([3]) ac- 
cording to 

i?(L,T) = J L 1 +"V / dr dr / G 2 (x,r;x',r'), (4) 
xx' Jo Jo 

the intersection of the curves R(Li,T) and R(Lj,T), 
shown in Fig.0 gives a size-dependent estimate Tl^l, {p) 
for the critical temperature T c (/j,) (26[. As L — > 00, the 
series of Tl^l^/j) converges to T c (/j) and one can ana- 
lyze it using corrections to scaling, to extract its limiting 
value ll|. Likewise, a linear fit of a size-dependent es- 
timate for the filling factor z/(L; /j) versus 1/L yields the 
thermodynamic limit filling factor v(fJ-). 

The next step is to repeat the procedure for a se- 
quence of fi values and extrapolate the resultant series 
of T c {v) towards v — > using the leading order form 
T c (v)/£p{i>) — T c /ef ~ const • v x l 3 . This functional form 
is expected from the analysis of the difference between 
the scattering T-matrices on the lattice and in the con- 
tinuum [m. 

Shown in Fig. [3] are the simulation results for the criti- 
cal temperature at filling factors ranging from 0.95 down 
to 0.06. We use system sizes up to 16 3 sites with up to 300 
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FIG. 3: The scaling of the lattice critical temperature with 
filling factor (circles). The errorbars are one standard devi- 
ation. The results of Ref. [23|, [24|] at quarter filling are also 
shown for a comparison. See the text for discussion. 



fermions. It is clearly seen that starting from v w 0.5 the 
expected v x l z scaling holds very well and the subleading 
corrections are negligible. On the other hand, close to 
half-filling, T c (y) is essentially constant, (see, e.g. |27|). 

Figure [3] shows a strong dependence of TJv) on v. 
This is in apparent contradiction with Ref. [161 ] which 
assumes no such dependence. This mi ght be due to 
different single-particle spectra : Ref. [161] employs a 
parabolic spectrum with a spherically symmetric cutoff, 
while we use a tight-binding spectrum over all of the 
Brillouin zone. Our preliminary tests show that ~ i/ 1 / 3 
corrections do depend on the specific choice of single- 
particle spectrum, and may even have different signs for 
different e/.. 

The critical temperature we derive from Fig. [3] is T c = 
0.152(7)ef- Various approximate schemes have in the 
past yielded T c to be either above [1, @, @| or below 0, @, 
G3 the BEC limit T BE c = 0.218£ F . Our results clearly 
show that it is below. 

Previous numerical results were also in disagreement 



on whether T c is higher or lower than Tbec : Ref- 1Z| 
quotes Tc/ef = 0.05, but the scattering length has 
not been determined precisely. Most probably, this re- 
sult corresponds to a deep BCS regime, where the crit- 
ical tem per ature is exponentially suppressed. Lee and 
Schafer [l8| claim an upper limit T c < O.liEp. This 
result is based on a study of the caloric curve of a uni- 
tary Fermi gas down to T/ep = 0.14 for filling factors 
down to v = 0.5. The caloric curve of Ref. 18j shows 
no signs of divergent heat capacity which would signal 
the phase transition. We find it not surprising since at 
quarter filling T c (y = 0.5)/e F ~ 0.05, see Fig. H The 
simulations of Ref. [lj| , which are also based on a caloric 
curve study, yield T c = 0.23(2)£p. What this otherwise 



excellent treatment lacks is an accurate finite-size and 
finite-density analysis of the MC data. 

The value of T c determined in this work cannot be 
directly compared to the experimental result T c — 
0. 27 (2) sp [3j for a number of reasons. First, there are 
strong indications that a presence of a trap significantly 
enhances the transition temperature, see e.g. Ref. [8j. 
Second, the data analysis of Ref. Q relies on a mean- 
field approximate theory for relating the empirical and 
actual temperature scales. In this regard, it would be 
extremely interesting to see to what extent the results of 
Ref. 3 would be affected if a different theoretical scheme 
is employed for thermometry. 
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FIG. 4: The temperature dependence of the energy per par- 
ticle (upper panel) and chemical potential (lower panel) of 
the unitary Fermi gas. Red circles are the MC results, black 
dotted lines and blue dashed lines correspond to the Boltz- 
mann and non- interacting Fermi gases, respectively, the dot- 
fashed lines are the asymptotic prediction of Ref. !29j (plus 
the first virial Fermi correction), black triangles are the path 
integral MC results of Ref. [lq |. and the purple stars denote 
the ground-state fixed-node MC results [la ]. 

At unitarity, the thermodynamic functions acquire a 
self-similar form [28l ]. For example, for the free energy 
one has: 



F(T,V,N) = f( F \T/E F )NE Fl 



(•5) 



where N is the number of particles, V is the volume, 
and f^ F \x) is a dimensionless function. Eq. ([5]) allows 
one to express all thermodynamic potentials in terms of 
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energy per particle = E/Nef and rescaled chemical 
potential = /i/ef- The latter quantities are directly 
measurable numerically. 

An analysis similar to the calculation of T c yields 

E/(Ns F ) = 0.31(1), (6) 
n/e F = 0.493(14). (7) 

For the pressure P and entropy S, one than has 

P/(ne F ) = 0.207(7), (8) 
S/N = 0.16(2), (9) 

which follows from (JT)) and exact relations PV = (2/3)E 
and S = (5E/3 - fiN)/NT. Eqs. CO)-© are for T = T c . 

Shown in Fig. |4] are our results for the dependence of 
the energy per particle and chemical potential on tem- 
perature. In the high-temperature simulations we use 
system sizes of up to 32 3 sites with up to 80 fermions. 

As can be seen in Fig. [4l our results for both en- 
ergy and chemical potential approach values close to the 
fixed-node MC values [H] as T -> 0. For T/sf_^ 0.5 
our results are not far from the curve of Ref. [la]. As 
T/ep — ► oo, both energy and chemical potential ap- 
proach the virial expansion [29} at high temperatures. 

In conclusion, we have performed a determinant di- 
agrammatic MC simulations of a unitary Fermi gas by 
means of diluting the attractive Hubbard model. In or- 
der to extract the continuum gas behaviour we carefully 
treat both finite-size and lattice corrections. We have de- 
termined the critical temperature T c /e F = 0.152(7), the 
values of the thermodynamic functions at criticality, and 
the overall shape of the thermodynamic potentials from 
zero- to high-temperature regimes. 
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Erratum: Critical Temperature and Thermodynamics of Attractive Fermions at 
Unitarity [Phys. Rev. Lett. 96, 160402 (2006)] 

Evgeni Burovski, Nikolay Prokof'ev, Boris Svistunov, and Matthias Troyer 



The value of entropy given in Eq. (9) of Ref. [1] has an incorrect error bar. The correct value is directly related to 
error bars for energy and chemical potential mentioned in Eqs. (6) and (7) using exact relation S/N = (5E/3N—/j,)/T. 
Equation (9) should read: S/N = 0.2 ± 0.2, i.e. we can not predict the correct value of entropy at the critical point. 
The loss of precision in the entropy estimate does not affect other results in [1]. However, it precludes quantitative 
discussion of adiabatic ramp experiments presented in Ref. [2] based on matching entropies of the non-interacting and 
unitary gases. In particular, the numerical value of (T/Tp) a given in Ref. [2] is meaningless. 

We thank W. Zwerger for pointing out the issue with the entropy calculation. 
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